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ABSTRACT 

A finite element program utilizing an axisymmetric eight degree 
of freedom isoparametric quadrilateral element capable of handling 
both isotropic and anisotropic (assuming similar properties in tension 
and compression) materials was developed. The problem of a thick 
walled open ended cylinder subjected to internal pressure was 
analysed and the finite element solutions were compared to the exact 
analytical solutions to show convergence to the correct answer and to 
give an indication of the order of magnitude of errors which might be 
expected in other problems to which no exact solutions are known. 
Several problems which could have significance in orthopedics and in 
prosthetic design were analysed, utilizing idealized bovine femur as 
a material, to give an indication of the potential value of finite 


element analysis in these areas. 
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CHAPTER I 
INTRODUCTION 


1.1 BROAD SCOPE OF THESIS 

Prosthetic devices are used in orthopedics to help heal some 
bone fractures as well as to assist mobility when the bone's load 
carrying capacity is diminished through disease. There are basically 
two components involved in prosthesis design - the element to be 
inserted into or attached to the bone and the bone itself. Design 
problems of this kind have, in the past, been solved by ad hoc 
methods. Only recently have systematic studies been performed to 
delineate the problems encountered in prosthesis design. Some of 
these problems are: mechanical and medical behaviour of bone-metal 
interface; response of live bone to various mechanical loadings; 
healing characteristics of bones when the fracture surface is under 
stress; response of prosthesis to various loads; and loads on 
prosthesis during walking or running. This list is far from 
exhaustive. The general problem of prostheses design is considered 
in this thesis. After a review of literature the aims of this 


work are stated after which the specific problems are formulated. 


1.2 REVIEW OF PERTINENT LITERATURE 
The need for accurate determination of stresses and strains 
in anisotropic bodies, particularly in bones, is increasing with 


the development of new and better prosthetic devices. 
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The earliest mention of the mechanical significance of bone 
form was made by Galileo (1638), but its strength and other physical 
properties were not extensively studied until the latter half of 
the nineteenth century. Since then stress and strain in bones have 
been studied by considering sections of bones, recording stress- 
strain phenomena in models of bones, analysing stress-strain phenomena 
in intact bones, and by mathematical analysis. A good review of these 
works may be found in Evans' [1] book Stress and Strain in Bones. 

Recently, finite element methods have been applied to the 
stress analysis of human intervertebral discs by H. Kraus et al [2]. 
The basic theory of the finite element method is developed by 
Zienkiewicz [3]. This, along with some reference to the fundamentals 
of anisotropic elasticity described by Hearman [4] and an isoparametric 
constant strain triangular element (F.E.C.S.T.) developed by Dr. D. 
Murray of the University of Alberta, was sufficient for the 
development of the axisymmetric isoparametric element used in this 
work. Zienkiewicz states [3], and Tong and Pian show [5] that this 
element formulation has the necessary requirements to converge to 
the correct solution. 

The fundamental equations of elasticity may be found in 
numerous texts - see, for example Wang [6], and have been applied 
to the isotropic open ended cylinder problem. The anisotropic 
solutions were developed from these equations without the 


assumption of isotropy. Bird et al [7] experimentally determined 
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the Young's moduli for bovine femur cortical bone in the radial, 
tangential and axial directions and gave an indication of the range 
of Young's moduli which may be expected in bovine femur. 
McElhaney [8], among other results for the dynamic response of bone, 
gave Poisson's Ratio of bovine femur for various rates of strain in 
the axial direction. Other works, including those of Dempster and 
Coleman [9 ], Evans et al [10] and [11], Amtmann [12] and Lugassy [13], 
gave an insight as to the complexity of structure and composition of 
bones and therefore of the complexity necessarily involved in the study 
of their mechanical behaviour. 

As Fung states [14]: 

"New experiments which cover a wide range of extensions in 
different directions are needed in order to extend our knowledge 
of the stress-strain relationship for biological tissues". 

It might also be added that methods of utilizing this 


knowledge, once obtained, must also be explored. 


1.3 AIMS OF THESIS 

The aims of this thesis were threefold. Firstly a finite 
element program was to be developed for determining displacements 
and stresses in isotropic and anisotropic axisymmetric, linearly 
elastic bodies under axisymmetric loadings. This program was then 
to be tested on a problem for which an exact solution exists to show 
convergence to the correct solution and to give an estimate of the 


errors which could be expected for problems for which no exact 
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solutions exist. 

Secondly, although quite generally applicable, this method 
was to be applied to analyse stresses in idealized long bones due 
to the insertion and presence of prostheses and hence to show this 
method as a useful aid for future prothesis design. 

Thirdly, further work necessary to give improved results and 
to obtain a better 'tool' with which to analyse stresses in bones 
(or other anisotropic bodies) and help design prosthetics was to be 


suggested. 


1.4 DEFINITION OF THE PROBLEMS TREATED 

Chapter II contains a description of the formulation of the 
finite element and the development of the stress-strain matrix [C]. 
In Chapter III the results and discussion for the test problem of 
thick walled open ended cylinders under internal pressure are 
presented. The assumptions used are stated and discussed in Chapter 
IV, and the results and discussion for bovine femur under various 
loadings are given. The loadings used are internal displacement 
due to the insertion of an interference fit rod, with and without 
a shearing traction applied by the rod, and axial force applied to 
bone with a flat and a hemispherical cap fixed on one end. 
Conclusions and suggestions for further research are given in 


Chapter V. 
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The basics of the finite element method are given in 
Appendix I. The governing equations and open ended cylinder solutions 
used for comparison in the test case are developed in Appendix II. 


The program used is described and given in Appendix III. 
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CHAPTER II 
GOVERNING EQUATIONS AND ELEMENT FORMULATION 


2.1 GOVERNING EQUATIONS OF TORSIONLESS AXISYMMETRIC ELASTICITY 
The following equations define the elastic axisymmetric problem: 


(i) Strain-displacement relations: 
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ore oy Sia Sr apy Fi ergs et ee 


r YZ r " 
ae aera, = 0 (2.22) 
oT Oo A 
YZ z rz 
Src i a 


(iii) The equilibrium equations in terms of displacements (Neuber- 


Papkovich) [21] have been expressed as: 
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(v) Hooke's Law for anisotropy with the principal directions of 


anisotropy coinciding with the r, z, and 6 directions: 
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Boundary conditions for the specific problems considered are given 
later. 
The general solutions for displacements of points of an 


elastic body in the axisymmetric problem have the form: 


u=4(v - 1)8, - S{o8, + 28, + By) 
(2.6) 
w= 4(v - 1)8, - S08, + 28, + 8) 
where Bos BS and Be” are harmonic functions, one of which may be 
set to zero without loss of generality, and » is harmonic also. 
Finding suitable harmonic functions which satisfy all 
boundary conditions of a particular problem is very difficult apart 
from some trivial cases, and therefore finding exact solutions is very 
difficult. 
Because of the difficulty in obtaining exact solutions and of 


the inflexibility of resolving the equations for each set of 
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boundary conditions to be considered, a finite element method was 
chosen to obtain solutions to the problems to be investigated. 

The exact solutions for the open ended cylinder under internal 
pressure are fairly easily worked out, and this is done in Appendix 
Dis 

2.2 FINITE ELEMENT FORMULATION 

The problem of stress distribution in bodies of revolution is 
essentially two-dimensional. By symmetry, the two components of 
displacements in any plane section of the body define completely the 
state of strain and, therefore, the state of stress. 

For the 8 degree of freedom quadrilateral element, with nodes 
i = 1, ..., 4 number in the anti-clockwise sense, a local 
coordinate system (&,n) was chosen (see Figure 1). 


The general internal displacements i and V are related to the 


nodal displacements u and v by: 
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where [NJ is the matrix of displacement coefficients. 


The general internal positions f* and Z are related to the nodal positions 


r and z by: 
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where [Nel is the matrix of geometric coefficients. For an 


isoparametric element such as this: 


{9}, = {lg = {91 
(2.9) 
or = IN,.] = ENg] = EN] 
For this element it can be shown that 
6, = TAC Sia] + ayer = 1, -.. 4 tice 10) 


where EO = € E. 
peat aati 
and E. and n; are values of &€ and n at node i. 
The strain vector defined below lists all non-zero strain 
components possible in an axi-symmetric deformation and defines them 


in terms of the displacements of a point. 
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Note: [Jo + ] is the conjugate transpose of [J]. 
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Jacobian and the comma indicates differentiation, 
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Theretore, substituting 2.03, 2.14, 2.16.and 2.17 into 2.11; 
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where [L] = [{o,,Hose}! - (s,.Ho.,)] 
. (€} = (81th 
where [B] is the element strain-displacement matrix. From the 
generalized Hooke's Law: 
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Expanding 2.20 and noting that for cylindrical coordinates: 


[ey .€5 29 2&4 Eo 2&6] * [e, rE 2E 9 2VonV29 Vg 


and [94405 203 104 05206] = [o, =) r2%q2T op rT 292 TRG 


and also noting that for axial symmetry: 


Yoal - Vest =T , = 0 (2522) 


The expanded form of Hooke's Law is then: 


(1/E,)o, ; (vi p/E,) op - (Vi9/E,)%9 
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where G = Gr = E,/2( + vee) 
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Ne ViglE, PVE N/E, 0 
0 0 0 1/G 
Therefore {6} = [C]{eE} (2.25) 


where [C], the constitutive matrix, can be found by inverting 
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For an isotropic material: 
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The elemental stiffness matrix, noting that the volume integral must 
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which may be written as 

[k]® = an | [8] [c][B]r DETLU]dedn (2.29) 
as drdz = DET[J]d&dn . 

As [B] and [J] depend on the coordinates, the Gaussian Quadrature for 


numerical integration 


n 
E HGHSF(E, on4) (2.30) 
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-] 
was used where for this case n = 2, H; = Hs = 1 and b, and n. take 
values +1/Y3. For n sampling points a polynomial of degree 2n - 1 
can be constructed and exactly integrated [3], therefore in this 
case the integration is approximate. 


The structural stiffness matrix [K] is obtained by the direct 


stiffness assembly of the elemental stiffness matrices: ; 


j.e.: [Ky] = r[K, 5° (2.31) 


where the contribution of each element to [K] must be evaluated 
individually (see reference [3].) 


Nodal displacements {43 are found by Gaussian elimination from 
{R} = [K]{4} (2.32) 
after the modifying of [K] and {R} for displacement boundary conditions. 
After the nodal displacements are found, stresses are computed 


for each element in turn from: 


fo} = [c] [Bb] . (2.33) 
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CHAPTER III 
TEST CASE 


In order to show that the program developed in Chapter I] 
gives reasonable results, and also to estimate the accuracy and 
convergence of the results, solutions were obtained by modeling open 
ended cylinders of various heights and internal to external radii 
ratios subjected to internal pressure and made of steel and bovine 
cortical bone. The properties of the isotropic steel cylinders used 
were E = 30,000 K.S.I. and v = .3. The properties of the bovine 
cortical bone are given in the discussion of assumptions in Chapter IV. 

These particular problems were chosen as test cases because 
exact solutions to them are available for the isotropic [6], as well 
as anisotropic behaviour (developed in the appendix). 

The finite element results for the principal stresses at the 
midplanes of these cylinders are plotted in nond imens ional] form in 
figures 2 to 19. These were all obtained using a 300 element grid 
(see Figure 1). 

3.1 TEST (CASE RESUIES 

Figures 2, 3, and 4 are qraphs of o./P vs r/b for isotropic 
steel cylinders of H/a = .8, 4, and 16 respectively, each having plots 
for cylinders of a/b = .25, .5, .67, and .91. The solid lines 
plotted are exact open end cylinder solutions corresponding to the 
finite element solutions at nodal points along the cylinder midplanes. 
The agreement between the exact and finite element solutions can be 


seen qualitatively to be good, the average error in radial stress: 
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Figure 1.b Three Hundred Element Grid. 
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at the nodal points being approximately .5% for cylinders with a/b 
from .25 to .67. The average error for a/b = .91 is higher: 7.88%, 
3.21% and 1.93% for H/a = .8, 4, and 16 respectively (see Table 1). 

Figures 5, 6, and 7 are similar to 2, 3, and 4, with the 
cylinder material being bone.. Here the average error is slightly 
larger than that in figures 2, 3, and 4, being approximately .6% for 
cylinders with a/b from .25 to .67. The average error for a/b = .91 
is 2.52%, 1.04%, and 1.32% for H/a = .8, 4, and 16 respectively, 
higher than that for other a/b ratios, as was the error in these 
cases for steel. 

Figures 8, 9, and 10 are graphs of p/P vs r/b for isotropic 
steel cylinders of H/a = .8, 4, and 16 respectively, the solid line 
plots again being exact open end cylinder solutions. Again the 
agreement between the finite element solutions and exact open end 
cylinder solutions can be seen to be good, the average error in hoop 


stress: 


being approximately 1.0% (see Table 2). Exceptions to this are for 
a/b = .91, where ey = 38.62%, 1.3% and 12.95% for H/a = .8, 4, and 
16 respectively and for a/b = .67 and H/a = .8 where Qy = 4.35%. 
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Table 1. Average Error in Midplane Radial Stress 


Table 2. Average Error in Midplane Hoop Stress 


1.03% 
1.88% 
4.35% 
38 .62% 
1.23% 


1.01% 


2.31% 
33.57% 


5a 


er 


B8e 0 
yan. 0 


wah. 0 


pee. 
928.0 
e289. 
gon.0 
oe 


$18.0 
met 


2eavd2 qooH anetgbtM nt sore si 


SISSe 
$°T.0 
efba, | 


or | 


my 


+ [ree 


: has 
ot bia 
We Brits 


ee. 


*saieputTA9 otdos4ostuy Joy sseaqs 


q/4 O ° if Q e 9 ° 
OY J : 3 6 
O-\ > S 
uy \y = 
Om wo x 
AY 
x 
¥ | 
4 
y 
4) 
\) 
Ay 
4) “ 
i \y 

4) 

i) 

4) uy 
43 
u 


Tetpey sue ldpth ° 


i 


SUOTYNTOS 4oexe 


iW 


- 


CG. 


BMS 


i) 


| 


ut 


SG eins 4g 


q/e O 
q/e O 
a/eV 
q/e O 


e/ 


of 


SAT Pinjete* 


acd 


ro 


£ 


ioe 


*saeputtsé9 otdozqostuy Jog ssezye Tetpey eueTdptyW °9 eansty 


oe Opt g° 9° 
~ = = : EEE - = ; 
YW WY a 
Vv 8 
iy \} &) 
) ‘. v 
Ly 
v v 
4) 
a 
(T) iv 
iy 
@ \ 
4) : \ 
rh 
4) 
4) WW 


% 


Y 


us 


(\) 


ps 


SUOTINTOS JOexeE 


T6* = a/e O 
49° = 4/e O 
Gc =q/fe Vv 
Se* = q/e O 
| ey 
> 


‘F.¥ D> 


€@xeGE 207 eT oxle 


24. 


*sioeputTAg oTdoujostuy Joy sserzqs [Tetpey eueTdpty °2 sans ty 


q/a Or g° 
ON, : 
‘VL 

iy Ve 
4) a 
ay . 
Y 
g &) 
4) 
rT 
iW 
4) 
1 
uN 
4 
4) 
ul 
y 
(} ) iW 


oe 


o) 


Vy. 


\" 


T° 


SUOT4NTOS JOexKE 


rm) 


(\) 


© 


q@ 


Ge 


Ow 


q/eO 


q/eo 


te 
\O 
i} 


Ge = gfe V 


Ww 
N 
iH] 


q/e oO 


OL = &/H 


v0 


4) 


OTL 


a/0= 


TIISOL 


a 


i 
th 
RA pe 
ou 
fh 


‘ 
ae 
” 
Le 
a 
os 


e <= 


Fou | 


opie 


26, 


g/t 


*sadeputTAg otdos,os—T Joy sseaqg dooy sue Tdptyl 


OL 


g 


T° 


SUOTYNTOS J0exe —— 


16". 


es 


i] 


°Q ouNnsty 


SO) 
q/e oO 
q/eV 
q/e Oo 


e/H 


G/T 


‘b 
ys) 
ee 
et 
i 

a 
Pe 


20). 


*szeputtAé9 otdorqosy— dog sseaqag dooy suetdptw °6 aanst gq 


Osi oh 9° ne Ge 


SUOTINTOS JOeXo —— 


T6° = 9/e QC) 


49° = 4/e Oo 


i) 
TT 
Q 
— 
oO 
<q 


OT 


> —— “exsop edpiproue 


ra 
— 


=@3.., 


QO \b =" 


oh. 


*sgeput{TA9 ofdotzosyT 40g sset4g dooy oueTdptn “OL eansty 


9° 


9° 


ine 


SsuOTyNTOS 


16" 
105 

G° 
GGe 


OL 


1oexe J— 


Ge 


a/e 


© 


Oo 
Vv 
Oo 


Gay 


eo 


ae! 
my 
wa 
& 
wo 
+3 


=i 
? 
t 


RTGDT aS yoo 


Fe" 


28. 


The reasons for these large errors are given in the following section. 

Figures 11, 12, and 13 are similar to 8, 9, and 10 with the 
cylinder material being bone. Again the magnitude of errors for bone 
are very close to those for the corresponding cases with steel 
cylinders (see Table 2). 

Figures 14, 15, and 16 are graphs of o,/p vs r/b for the same 
sets of steel cylinder geometry mentioned in reference to figures 
2, 3, and 4. For open ended cylinders theory predicts that Oo, should 
be zero everywhere for loading by internal pressure only, therefore 
these graphs also indicate the error in o,/p. The errors in axial 
stress can be seen to be generally small when compared with p, and 
are largest close to the inner cylinder radii. As was the case with 
errors in radial and hoop stresses, the errors in axial stress are 
larger for cylinders with a/b = .91 than for other cylinders. 

Figures 17, 18, and 19 are similar to figures 14, 15, and 
16 respectively, with the cylinder material being bone. The errors in 
axial stress in these cases are of the same order of magnitude as 
those for corresponding cases for steel cylinders except near the 
inner radii where the errors, although still eialieare Significantly 
larger than those in the corresponding isotropic cases. 

For torsionless axisymmetry the shear stresses Ty, gand TO9 
are necessarily zero, and for the test cases Ty should also be zero 
everywhere in the cylinders [6]. In all cases the errors in 
nondimensional shear stresses Taal P predicted by the finite element 
program were very small (several orders of magnitude smaller than the 


corresponding errors in o/p) and were therefore not plotted. 
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Figures 20, 21, and 22 are graphs of -0,./P VS Z/Z9s 
g,, /P vs z/Zp and o,/p vs z/Z respectively for a steel cylinder of 
H/a = 4 and a/b = .5 at r/b = .5, .75 and 1.0; 7.e.the inside, middle, 
and outer radii of the cylinder. These results were obtained using 
the 300 element grid of Figure 1. The solid lines are again open end 
cylinder solutions. These plots indicate that there are increased 
errors in the finite element results at the ends of the cylinders, and 
give a rough idea of the magnitudes of this increase for this grid. 

Figures 23, 24, and 25 are similar to 20, 21, and 22 with 
the cylinder material being bone. These also show the presence and 
magnitude of the increased errors at the cylinder ends. 

Figures 26 and 27 are graphs of the log of the absolute 
percentage error in - o/P and o,/P (compared to the exact solutions) 
versus the log of the number of elements in the radial direction 
through the cylinder wall. A cylinder of H/a = 4 and a/b = 5 was 
used for this illustration of the convergence of these stresses to the 
correct answer [5]. Values of these errors were plotted for both 
isotropic and anisotropic (bovine bone) cylinders at r = a and 
r = (b-a)/4 at the mid-height of the cylinders. The errors in these 
stresses can be seen, in general, to decrease with an increase in the 
number of radial elements. This approximately follows a relationship of 
the form e= A/n* where n = number of radial elements. A and x are 
constants which vary with position and mechanical properties of the 


cylinder. 
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Figure 26. Radial Stress Convergence. 
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Figure 28 is a plot of absolute error in o,/P versus the log 
of the number of radial elements. Absolute error was used here instead 
of percentage error because the true axial stress in this test case 
should be zero, which makes a percentage error undefined. This figure 
indicates that the axial stress is converging to the correct solution 
in a similar manner as the radial and tangential stresses. 

Cylinders under internal pressure, modeled with varying 
numbers of elements in the axial direction only, indicated that for 
this loading the accuracy of results did not depend on the number of 
axial elements; (i.e. there was no convergence of any stresses for an 
increased number of axial elements). This was expected because ideally 
stresses should change very little in the axial direction (at a 
certain radius) under this type of loading. 

3.2 DISCUSSION OF TEST CASE RESULTS 

In general it would be expected that for the test cases used 
here, the accuracy of the results would increase as the a/b ratio 
approached unity. This is because as a/b approaches unity the radial 
and tangential stresses more closely resemble straight lines when 
plotted against the radius; the axial stresses are constant, and 
therefore the linear strain finite element used here can more accurately 
model the strain and stress distribution [3]. 

It can be seen from Tables 1 and 2 and Figures 8, 9, and 10 
that this is true to a certain extent. The error decreases for a while 


in all cases as a/b increases, but then increases. This is mainly due 
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Figure 28. Axial Stress Convergence. 
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to the number of elements used becoming too large, (i.e: for a completely 
linear strain or stress distribution one element should give very 
accurate results) increasing the numerical errors incurred. For the 
case Of a/b = .91 much more accurate results could be obtained using 
fewer elements. The 300 element grid was used throughout the test cases 
(except for convergence tests) because the problems discussed in 

Chapter IV were solved using this grid or grids similar to the 300 
element. grid. For these problems a/b was equal to .5, which is in the 
region of good accuracy for this grid. 

Some of this error increase could also be due to the increased 
"aspect ratio"[15] of the elements as a/b approaches unity. This is 
numerical error due to the lengths of the sides of an element being 
greatly different. As this must also account for most of the 
variation of error with cylinder height seen in the tables and figures 
mentioned above, it can be seen that the errors caused by it can become 
quite large, especially as a/b approaches unity, and particularly for 
the hoop stress. (See tables 1 and 2). 

Larger than average errors were found to occur near the 
boundaries of the cylinders, and in particular near the inner radius. 
This is probably partly due to the fact that loadings which should be 
distributed evenly over the inner surface are assumed concentrated on 
external nodes, and also that stresses and strains are changing most 
rapidly near the inner radius [3]. These errors were lessened by 


using more elements near the inner radius and at the cylinder ends. 
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The errors found for the anisotropic (bone) cylinders were 
generally slightly larger than the corresponding errors for isotropic 
cylinders at low a/b ratios and lower at high a/b ratios, however for the 
a/b ratio of the bovine bone model used in the problems of Chapter IV 
for this grid, the magnitudes of the average errors in stress are about 
the same. 


It should be noted that the'average errors' used here are defined 


as 


where n is the number of nodal points along the mid-height of the 
cylinder and e. is the absolute percentage error in stress at that node 
for radial and tangential stress and the absolute error for axial stress. 
This gives a weighted average as there are more nodal points near the 
inner radius. As the largest errors occur here, this ‘weighted average' 
is larger than true average would be. 

In Figures 26, 27, and 28, the errors in stress were seen to 
decrease as n increased, then increase as n further increased at 
r = (b-a)/2. The increase in error here was probably due to increased 
numerical error, ie: too many elements for optimum accuracy. At this 
radius the radial and tangential stresses change at a slower rate than 
near the inner radius, therefore fewer elements are required here to 
give the same accuracy as near the inner radius. 

The errors in the finite element results could in some cases 


be lessened by further subdivision of the elements. The limit of the 
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effectiveness of this subdivision is the point where roundoff errors 
become the same order of magnitude as the errors caused by the 
finite element size. After this point, greater accuracy can be 
achieved by using double precision and then possibly still further 
reduction of element size. Generally, for most applications 
sufficient accuracy can be obtained with single precision and a 
relatively coarse grid size. 

It should be noted that the large errors in hoop stress (e,) 
are in most cases due to the definition cf this error. By 
comparing the absolute error to the range of 0, across the 
cylinder wall, errors very small in comparison with the actual 
value of 0, result in large errors as the range is much smaller 


than the value of 0,: 
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CHAPTER IV 


RESULTS AND DISCUSSION OF THE FINITE ELEMENT 
ANALYSIS FOR BONE UNDER VARIOUS LOADING CONDITIONS 


4.1 ASSUMPTIONS 

In the problems analysed the following assumptions were made: 

1. The bone material is homogeneous. 

2. The portion of the bovine femur used is axisymmetric. 

3. The material has constant elastic physical properties. 

4. The material has similar physical properties in tension and 

compression. 

5. Rods and prosthetic devices are rigid. 

6. Strains are small compared to unity. 

The assumption of homogeneity of bone is incorrect in that bone 
is not homogeneous in the microscopic sense. Bone is composed mainly 
of light, dark, and intermediate osteons which are made up of collagen 
fibers, interfibrillar calcified substance, and pores (Haversian 
systems and Volkmann's canals). [12] and [10], The modulus of 
elasticity and mechanical properties of bone are dependent upon collagen 
fiber orientation and distribution, and on the composition and 
distribution of minerals in the interfibrillar calcified substance. 
Fortunately, however, these components are small enough that for 
engineering purposes, in most cases, a segment of bone when viewed 
macroscopically is virtually homogeneous in terms of its mechanical 
properties. In the macroscopic sense, although as just stated, a 


segment of bone may be considered homogeneous, different parts or 
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orientations in a bone, different bones in a particular animal (or 
human), and similar bones (and areas in bones) in different animals 
of the same species may differ substantially in mechanical properties 
and physical makeup, and the same part of the same bone may change 
with time. When considering one portion (say the middle third) of a 
long bone (cortical bone), however, the assumption of homogeneity in 
the macroscopic sense is fairly well justified in most cases [11]. 

The assumption of axisymmetry of a portion (middle third) of 
a typical long bone is obviously an inaccurate one, however, the 
assumption of any "typical" cross-section would be inaccurate, as 
the true cross-section of a particular bone will vary from place to 
place in that bone, and the cross-section of similar bones in different 
animals will vary also. Usually before the insertion of a prosthesis, 
the intramedullary canal (interior of the bone) is reamed out to a 
near perfect cylindrical bore, but the outside of the bone is 
untouched. This makes the assumption of axisymmetry (particularly 
near the interior of the bone) a little better. Despite the errors 
necessarily induced by this assumption it was felt that useful 
information could still be obtained from an axisymmetric bone model, 
in that results for the model would be of the same order of magnitude 
as those for a real bone particularly near the more axisymmetric 
interior of the bone. 

Constant elastic physical properties (Young's moduli and 
Poisson's ratios) were assumed for the bone models. As bone is 


viscoelastic, this is not absolutely true, the physical properties 
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actually being dependent on both strain, strain rate and previous 
stress-strain history [7],[13], however in general the assumption is 
acceptable as the stress-strain curve for bone is nearly linear up to 
failure [17]. The assumption of perfect elasticity itself is not 
quite correct, but in the range of stresses with which the test problems 
deal, it is probably close. The elastic constants used for the test 
problems were E. = 670 ksi, EY = 1500 ksi, EL = 2500 ksi, Ving = a = 
v_, = .3. These were taken from references [7] and [8], the moduli of 
elasticity being estimates for low loading rates of bovine bone from 
reference [7] and Poisson's ratio was taken from reference [8] (Vig : 
on and Vr were assumed to be the same as Vag aS no published data 
could be found giving their values. 

The assumption of similar physical properties of bone in tension 
and compression is probably fairly good. Comparison of moduli of 
elasticity in tension [11] with those in compression [7] for 
longitudinal cortical bone samples indicates that this assumption is 
good for this particular property. The values of E, in tension and 
compression are the same within the limits which could be expected 
due to the differences in different bones. It was assumed that this 
was true for the other properties as well, as no data could be found 
to make further comparisons. 

The assumption of rigid rods, wedges and prosthetic devices is 
justified in view of the fact that the modulus of elasticity for 
steel (of which most prosthetics are composed) is more than an order 


of magnitude greater than that of the largest modulus of elasticity 
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Infinitesimal strains were assumed for the finite element model, 
and therefore significant errors would arise if strains were allowed to 
become too large (greater than .005"/in). If results are converted 
from the non-dimensional form presented here it should be ascertained 
that the strains for a particular case are not beyond the limit 
mentioned above. 

4.2 RESULTS AND DISCUSSION 

To illustrate the usefulness of the finite element method in 
general and of this program in particular to the analysis of 
anisotropic axisymmetric bodies, specifically idealized long bones, 
two test problems were analysed: 

(1) Idealized long bones strained by the insertion of 
interference fit rods with and without shearing traction 
applied. 

(2) An idealized long bone under axial loading with flat and 
hemispherically shaped prosthetic caps. 

The first problem, that of interference fits of rods into bones 
has significance in orthopedic intramedullary nailing techniques. 
Here a rod is inserted in the reamed intramedullary canal of a long 
bone (usually broken) to increase its rigidity during healing. At 
the present ‘state of the art' the canal is reamed to a slightly 
larger diameter than the rod to be inserted and the rod is inserted from 
one of the bone ends and glued in place. It is possible, that in 
certain cases the use of an interference fit rod of shorter length 


and inserted into the bone halves through the break area could 


i i / Val rY iy ’ 


haSVaVAND S716. se tiiels ve “at\200. ia pike 
bantedyedes ed “pT borle aT oie bednbeong mio f si | 


ttme | ant bioysd Jon aie! ‘meen 2 astust 6g Cl Ot! 


| | - voreaiete/ dim an 3 ft ; 
At Bote $name ts otintt ont Me deontwiot tt ote bpoiiabiel ‘ . 
| 70 eiayions sto nfuotteg mt mengond Sry Ye. i 
<z9n09 poo! bashtgsbh eisai Tig9q2 . 2otbod. tSanmyetss |5t qo 
 bseutond on motdonq $29: 
to notdseent ‘ond on bentende eanod vier Cie tt ipa 
nottoe4 paisa vont iw bas tw abo a i 


uh ia =e Ai es 
\ ia : D ihe ¢ 
( 


Aye 
oe 


Tey 
i 


56. 
achieve greater rigidity and cause less trauma to the rest of the bone 
(in particular to the bone's blood supply). This problem may also have 
Significance in the insertion of some prosthetic devices (e.g. knee 
joints etc.) into bone, where a rod is pushed into the intramedullary 
canal to help secure the prosthesis. Here some interference fit may 
help stabilize the prosthesis and avoid localized contact areas 
between the prosthesis and the bone. [14] 

Figure 29 shows the shape of an idealized long bone of a/b = .5 
before and after the insertion of a rod with an interference fit of 
I/a = .005, where I is the difference between the rod radius and the 
inner bone radius. 

The boundary conditions for this problem are; u/a = .005 from 


z/z. = 1/3 to 1 on the inner radius, where u_ is the radial 


0 
displacement; v/a = 0 at 2/Z, = -]1, where v is the axial displacement. 
All other boundaries are free. 

The outline of the top half of the bone before insertion of the 
rod is drawn in heavy lines, while the outline of the bone after 
insertion of the rod is drawn in light lines. The non-dimensional 
displacements are exaggerated by a factor of 10° (i.e. an actual 
displacement of .01/a would be displaced by 1 inch on the drawing 
in both the r and z directions). 

This figure shows the radial displacement of the outer radius 
of the cylinder as nearly constant near the top of the cyliner, and 
approaching zero below the lowest part of the rod in a gentle "S" 


shaped curve. The radial displacement of the inner radius of the bone 
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Figure 29. Bovine Femur Shape After Rod Insertion. 
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corresponds to the outside of the rod to where the rod ends and 
approaches zero below the rod in a sharp curve. This curve is 
indicated as a dashed line as the grid used here was too coarse to 
show the "S" nature of the curve near the edge of the rod. The 
actual curve must be an "S" curve to prevent a discontinuity in the 
first derivative of the displacement. The cylinder is also shown to 
shorten in height somewhat due to the insertion of the rod. 

Figure 30 shows the shape of the bone of figure 29 before 
(dark lines) and during (light lines) the insertion of the interference 
fit rod. The boundary conditions here were the same as those for 
figure 29 with the exception that at r = a, from z/Z5 = 1/3 tola 
shearing traction was applied. The magnitude of the shear was 
determined assuming a coefficient of sliding friction of .25 between 
the rod and bone, and a normal force due to 0. found from internal 
displacement alone. Only this one iteration was taken as the coef- 
ficient of friction is assumed anyway. The non-dimensional radial 
displacements are again exaggerated by a factor of 10°. The radial 
displacements are similar in this case to those shown by Figure 29 but 
due to the shearing traction the axial displacements here are much 
greater. The inside upper edge of the cylinder has been 'rounded' due 
to the shearing traction and the top of the cylinder slopes down towards 
the centre of the cylinder. The axial displacements are only exaggerated 
by 10. 
The non-dimensional radial and hoop stresses (0,/0,) for the 
bone of Figure 29 pathoug shearing traction (solid lines) and for the 
bone of Figure 30 with shearing tractions applied (dashed lines) are 


plotted in Figure 3] versus r/b where 0, is the radial stress at 
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Figure 30. Bovine Femur Shape During Rod Insertion. 
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Radial and Hoop Stresses During 
and After Rod Insertion. 
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r = a and z/Z, = .5 for I/a = .005 and no shearing traction applied. 
These are plotted for values of z/Z, =702833), and: <5. 

Figure 32 consists of plots of non-dimensional axial and 
Shearing stress versus r/b for the problems mentioned above. In this 
case Tt, is the shearing traction applied at r = a and 2/2 = 5 for 
I/a = .005 and up = .25. These are also plotted for values of 
z/Z, =O, F so vanG:. 6O: 


Figure 33 is a series of plots of o versus I/a for 


emax/ Teult 
idealized long bones of various a/b ratios, where Soult is the ultimate 
tensile strength of bovine bone in the hoop direction. No shearing 


traction is applied here. The maximum hoop stress ( ) was found 


° emax 
in each case to occur near the tip of the rod, i.e. at z/Z = 1/3. 

The maximum hoop stress was plotted because it was felt that for this 
particular problem it is an important (and perhaps the most important) 
criterion for failure. To predict adequately the failure of bone 

a failure theory for anisotropic bodies probably taking into account 
elastic deformation due to hydrostatic pressure or tension, would 

have to be developed: These plots indicate that the maximum hoop 
stress varies linearly with the ratio of interference fit to internal 
radius of the bone. 

Figure 34 is a series of plots similar to those in Figure 33 
with the exception that here a shearing traction is applied to the 
inside radius of the bone as if the rod were being pushed into the 
bone. The magnitude of the shear was determined assuming a coefficient 


of sliding friction of .25 between the rod and bone. Again the 
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O Axial Stress at 2/ Zo = 0 
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Figure 32. Axial and Shear Stresses During 
and After Rod Insertion. 
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Infinitesimal Strains 
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Interference Fit After Rod Insertion. 
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Figure 3. Hoop Stress vs Interference Fit During Rod Insertion. 
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maximum hoop stress was found to occur near the tip of the rod in 

each case and the maximum hoop stress varied linearly with the ratio 
of interference fit to internal bone radius. The magnitude of the 
maximum hoop stress under these conditions (shear applied) was found 

to be less than the magnitude for the corresponding case without shear. 

Such plots as those which occur in Figures 33 and 34 could 
possibly be drawn up for various bones, sizes of bones, similar bones 
in different age groups etc. and could be beneficial in orthopedics; 
providing easy to use, quickly accessible estimates of maximum bone 
Stresses. 

The second problem is that of a bovine bone with a prosthetic 
rigid cap fixed on one end (as if it were glued). An axial force is 
applied causing an axial stress at the interface of the bone and cap. 
This is similar to the situation found at the bone prosthesis interface 
of an artificial knee joint when weight is placed on the leg as in 
standing. In the first case a flat rigid cap was analysed and in the 
second case a cap hemispherical across the bone wall was analyséd. * 

The boundary conditions for these cases are: 

at z/Z, = 0, oS P, 
at the bone-cap interface u = v = 0 
all other boundaries are free. 

Figure 35 contains plots of o/P, and op/P, versus r/b at the 
bone-cap interface and at 2/2. = .9 for the flat cap (solid lines) 
and the hemfepherteawwean (dashed lines). Both the radial and hoop 
stresses approached zero as z/Z, approached zero. 


* The toroidal cap is hereafter called a hemispherical cap. 
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Figure 35. Radial and Hoop Stress for Capped Bone. 
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Figure 36 contains plots of ash. and tee hD, versus r/b at the 
bone-cap interface, at z/Z, =.9,. — Again the.solid 
lines indicate the flat cap and the dashed lines indicate the 
hemispherical cap. This figure indicates that there are axial stress 
concentrations at the inner and outer radii of the bone for both caps, 
but that the average axial stress at the interface is less for the 
hemispherical cap than for the flat cap. This is due to some of the 
axial force being taken up in the form of shear stress along the 
interface in the case with the hemispherical cap. 

Figure 35 indicates that the average radial and hoop stresses 
at the interface are smaller for the hemispherical cap as well as 
the axial stress. 

Plots such as these could be helpful in the design of 
prostheses, pointing out areas of high stress and giving an 
indication of the best prosthetic shapes to use to eliminate stress 


concentrations or to lower average stresses in problem areas. 
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Figure 36. Axial and Shear Stress for Capped Bone. 
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CHAPTER Y 


CONCLUSIONS AND FUTURE RESEARCH 

».1 Conclusions 

A finite element program capable of modelling axisymmetric 
linearly elastic bodies under axisymmetric loadings was developed. 
This program was shown to give reasonably accurate results for both 
isotropic and anisotropic cylinders subjected to internal pressure. 
The program was applied to the solution of two problems, the results 
for which are given in Chapter IV. The accuracy of these results 
should in general be as good as the accuracy shown for the 
anisotropic cylinders under internal pressure as similar grids, 
cylinder geometry and elastic properties were used throughout. Larger 
errors than normal may be expected to occur at points where stresses 
are changing rapidly (e.g. at the tip of the rod in problem 1) and 
where singularities exist (e.g. at r = a or b and Lye ear On 
problem 2.) The results in these instances should be taken with caution. 

The results obtained in Chapter IV show that within the limitations 
imposed by assumptions, this method could be useful as an aid for 
future prosthesis design and in the analysis of stresses and strains 
in long bones. As most real bones are not adequately modeled under 
the assumptions used here, more research is needed to eliminate the 


need for as many assumptions as possible. 
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2.2 Future Research 

To model adequately the geometry of some parts of Jong bones 
and of many other bones, a three dimensional anisotropic element 
(e.g. a curvilinear hexahedron) should be developed. As well as 
adequate geometric representation, such an element would be able to 
model non-axisymmetric loading conditions and could to a certain 
extent account for regional differences in the material properties 
of bone. This is because each element may be given different 
properties. 

More information on the material properties of different bones 
and positions in bones is needed. In particular, estimates of Young's 
moduli and Poisson's ratios for the radial, tangential and axial 
directions in long bones are required to predict accurately stresses 
and strains in these bones. 

Present commonly used failure theories (Tresca and Von Mises) 
[18] do not take into account hydrostatic pressure which may produce 
distortion in anisotropic materials. A failure theory for anisotropic 
materials should be developed to aid in the design of prosthesis which 


will not cause bone failure. 
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APPENDIX I 


BASICS OF FINITE ELEMENT DISPLACEMENT METHOD [3] 

Conventional engineering structures can be visualized as an 
assemblage of structural elements interconnected at a discrete number 
of nodal points. If the force-displacement relationships for the 
individual elements are known, it is possible to derive the properties 
and study the behaviour of the assembled structure. 

In an elastic continuum there are an infinite number of 
interconnection points and here lies the biggest difficulty in 
obtaining a numerical solution (e.g. for stresses or strain) for the 
continuum. The concept of finite elements attempts to overcome this 
difficulty by assuming the real continuum to be divided into elements 
interconnected at only a finite number of nodal points at which some 
fictitious forces, representative of the distributed stresses actually 
acting on the element boundaries, are supposed to be introduced. This 
idealization reduces the problem to that of a conventional structural 
type which may be treated numerically. 

Using the displacement approach, this may be done in the 
following manner: 

a) The continuum is separated by imaginary lines on surfaces 

into a number of 'finite elements’. 

b) The elements are assumed to be connected at a discrete 

number of nodal points situated on their boundaries. 
The displacements of these nodal points, {6}, are the 


basic unknown parameters of the problem. 
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A function (or functions), {oh is chosen to define 
uniquely the state of displacement within each 'finite 
element' in terms of its nodal displacements. 


e.g. for the case {6} = {} 


oom 

£2 

Ww 
iT] 


(o}! {u} and 


So 

<2 

ees 
iT] 


{oH Wy) 
where {ti} and {V} are general internal displacements 
and {u} and {v} are nodal displacements. 
This function must impose compatible deformations on the 
structure,1.e. always satisfying geometry and boundary 


conditions. 


Note: For an isoparametric element toh, = {o}, = {o} 


d) 


where to}, is a function or functions giving the 
relation between coordinate locations in the deformed 
element. e.g: {if} = Cohein} 
| {2} {o}4tz} 
where {*} and {2} are the coordinate locations in the 
deformed element. In the above equations 
{r} and {z} are the nodal coordinates of the deformed 
element. 
The displacement functions now define uniquely the state 
of strain within an element in terms of the nodal 


displacements . . 
{é} = [B] {6} (A1.1) 
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with any initial strains and the elastic properties of 
the material will define the state of stress throughout 
the element and hence also on its boundaries. 

e) A system of forces concentrated at the nodes and 
equilibrating the boundary stresses and any distributed 
loads is determined, resulting in a stiffness relation- 
ship for the structure. 

To make the nodal forces statically equivalent to the actual 
boundary stresses and distributed loads, the simplest procedure is to 
impose an arbitrary nodal displacement and to equate the external and 
internal work done by the various forces and stresses during that 
displacement. This is the principle of virtual work. 

Let such a virtual displacement be d{6}© at the nodes. This 


results in displacements and strains within the element equal to 


d{f}= [N]d{o}® (A1.2) 
and d{é} = [B]d{o}® (A1.3) 
respectively as: Se A 
(£}= ENJC6}® = ENjy Njs Nye Df (A1.4) 
m 


= vector of displacements at any point 


within the element. 


and {é} =/e.,} = [B]{6}°= vector of nodal strains. (A1.5) 
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Lore 
where LN oan o! (Al .6) 
and superscript e designates elemental properties. 
The work done by the nodal forces is equal to the sum of 
the products of the individual force components and corresponding 


displacements, 


fe (enone) baat (A1.7) 
a 
; 
where {F}E = F 


are nodal forces statically equivalent to the boundary stresses and 
distributed loads on the element. Similarly the internal work per 


unit volume done by the stresses and distributed forces is: 


d{é}!{o} - dis}! {p} (Al .8) 
where {o} is the stress vector = [C]{€} (no initial strains or stresses) 
{p} is the vector of distributed loads acting on a unit volume 
of material within the element with directions corresponding to 


those of {6} at that point, 


or (d(s}°) "([B] {0} - [N]"{p}) | (A1.9) 
Equating the external work with the internal work obtained by integrating 
over the volume of the element: 
Q,T--1e@ _ e,T T qT (Al .10) 
(d{6}®) "eF3® = (dto}®) '(f[B] {o}d(VOL) - f [IN] {p}d(VOL)) 
As this relation is valid for any value of the virtual displacement, 


the equality of the multipliers must exist. 


(F}® = ({[B]"[C][B]d(vOL)){6}° - SIN] '{p}4(VOL) (Al.11) 
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This relation can be recognized as one typical of characteristics of 


any structural element in the form: 


{F}° 


[K]°{6} + {FY (A1.12) 
Comparing the two: 


[k]® = f[B]'[C][B]d(VOL) = the element stiffness (A1.13) 
matrix 


and nodal forces due to distributed loads are 


(FES = - JIN] "{o}d(VoL) . (A1.14) 
In general, external concentrated forces may exist at the nodes and 
the matrix | | 
Ry 
tr} =| 2 
R 
n 


will be added to the consideration of equilibrium at the nodes. 

If at the boundary, displacements are specified, no special 
problem arises, however, if the boundary is subject to a distributed 
loading, say {g} per unit area, a loading term on the nodes of the 
element which has a boundary face will have to be added. By the 


virtual work consideration, this will simply result in 


(Fie = -f [NI] '{ghd(area) (Al .15) 
An integration of this type is seldom explicitly carried 
out. Often by ‘physical intuition’ the analyst will consider the 
boundary loading to be represented simply by concentrated loads acting 


on the boundary nodes and calculate these by direct static procedures. 


(SF. 7A) 
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Once the nodal displacements have been determined by 
solution of the overall ‘structural’ type equations, the stresses at 


any point in the element can be found from: 


for SCC iLer cor (A1.16) 
This development for an element may be applied directly to 


the whole continuum, that is let: 


{f} = (N]{s} (Al.17) 
in which {6} list all the nodal points and N = Ne when the point 
concerned is within a particular element e and 7 is a point associated 
with that element. If point i does not occur within the element, 
N. = 0. Matrix [B] will follow a similar definition and the virtual 
work principle will be now applied to the whole structure. Interelement 


forces no longer need be considered and external virtual work during any 


virtual displacement of all nodes d{é} becomes 


d{S}!{R} (A1.18) 


while the internal virtual work is 
[ace) "ota voL) i [ac"ep)a(vou a facea"tord (avery (A1.19) 
Vv V 


where the integrals are now taken over the whole region. 
On substitution of 
d{f} = [N]d{é} and d{é} = [B]d{6} 
together with {fo} = [C]{é}. 
we have immediately on equating internal and external virtual work 
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and typical terms of the stiffness matrix become 

[K; 5] = §[8],(C][B,J4(voL) (Al .21) 
with the integral being taken over the whole region. 

Considering, however, the relationship between [B], and [B]. 
it is obvious that 


[K; 3] = Lk; 3° (A122) 


The same can be shown to be true of the various force components. 
The principle of virtual displacements ensures satisfaction 
of equilibrium conditions within the limits prescribed by the assumed 
displacement pattern. If the number of parameters of {6} which 
describes the displacement increases without limit then ever closer 
approximation of all equilibrium conditions can be ensured. 
The virtual work principle can be restated in a different 
form. Equating (18) and (19) : 
id{é}"{o}d(VOL) - [d{o}"R + [eA ¢p}a(vou) + [aca "ig}a(area)] =0 
} 


(Al .23) 
The first term of (23) is the variation of strain energy of the 


structure while the second is the variation of the potential energy 
of the external 10ads. 

Equation (23) means that for equilibrium to be ensured the 
total potential energy must be stationary for variations of admissible 
displacements. The finite element equations developed previously are 
the statements of this variation with respect to displacements 


constrained to a finite number of parameters {6}. 
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It can be shown that in elastic situations the total 
potential energy is not only stationary but is a minimum [16], Thus 
the finite element process seeks such a minimum within the constraint 


of an assumed displacement pattern. 
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APPENDIX II 
ELASTIC OPEN ENDED CYLINDER SOLUTIONS 
Considering cross-sections that are far enough away from the 
cylinder ends such that planes before cylinder distortion remain 
planes after distortion, stresses and strains are functions of the 


radial coordinate only. Therefore equation 2.4 becomes: 


deg 
Raat peace. ~ 0 (A2.1) 
As the cylinder is open ended, cuales 0 everywhere. 
Therefore, from equations 2.5 
Ee = C159 + C60), (A2.2) 
ey epg? 3%y 
where C, = T/Eps Co =-V.g/EVs C, = T/E,, 
al et Sa. 10), 
dr T dr a. av 
Substituting into equation (A2.1) gives: 
do do 
peu is S 2 2 ; 
r(C, oe Cy a ) + C15, + C40, Cop C30,, Qos (A2.3) 
Choosing 
epliayeI 
Oe ea ey 
2 
Oy = oye y! (A2.4) 
or 
do 
reas: ] 
die & ne. . Ta 
do 
0 _ i 
dr ey 
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the equilibrium equations (2.2) are satisfied identically. 


Substituting (A2.4) into (A2.3) and collecting terms gives: 


yl tay! = te ‘ (A2.5) 


Petco" In-r | (A2.6) 


C 
oy -ey=0 (A2.7) 


which is an ordinary differential equation with constant coefficients 


to which the solution is: 


We Gate ge 4 Aye’ as | (A2.8) 
where Cy = C2/C, and A, and A, are constants 


or, using (A2.6) 
y=A, en me Aor WC, (A2.9) 


For an open ended cylinder under internal pressure only the boundary 


conditions are: 


at r=a, o, = =y(a) = -p 


(A2.10) 
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where p is the internal pressure. Substituting these into (A2.9) 


gives: 
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(A2 


(A2 


17) 


elz) 


13) 


“ap = A, aya 4 eae) 
O=A, pia, foo Ca 
Multiplying (A2.12) by aay “4 and subtracting this from (A2.11) 
yields: 
A, =(-ap b"~4)/p 
where p= arg pea _ yvCag-vCq 


Substituting (A2.13) into (A2.12) and solving for A, gives: 


vCq 


A )yDbe 


9 =(ap b 
Substituting (A2.13) and (A2.14) into (A2.9): 


ap 


y= BB (EE. MCE ) 


Therefore from (A2.14) and (A2.15) 


8 (gH VCH _ 4VCH, VA) 


Oy 2 sete, -vC4 pVe4vCa 


From equations (A2.2) Cy = E/E, 


For isotropy Cy = 1 and (A2.16) and (A2.17) reduce to the tamé 


equations: 


2@i2 
a. = ap (Bo) 
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APPENDIX ITI 
FINITE ELEMENT PROGRAM 
This program performs a finite element analysis for linear, 
elastic axisymmetric problems using an eight degree of freedom 
isoparametric quadrilateral element. It outputs displacements, 
element stresses, element principal stresses, maximum principal 
stresses , nodal stresses, non-dimensional nodal stresses, and nodal 
corrections due to weight of element material. The flow chart, 
showing all subroutines, is presented in Figure 37. 

3.1 Explanation of Subroutines 

MAIN - Calls other subroutines in order. 

READIN - Reads and prints input data. Generates intermediate 
nodes and elements. Calculates half-band width and 
number of equations. 

ASTIF - Forms the element stress-strain matrix if one 
material only is used. Assembles the element 
stiffness matrices into STIF. Assembles the applied 

Toad vector. 

ELSTIF - Forms the element stress strain matrix if more than 
one material is used. Forms the element stiffness 
matrix for each element in turn. 

BETA - Forms the element strain-displacement matrix and 
element stress matrix for each element in turn. 

MODIFY - Modifies the stiffness matrix and applied load 


vector for displacement boundary conditions. 


| 8 


eeeAntd | 
— mobsoyy 40 pit Toten ts ens a . 
- atdemsosiqzt re adugdud a, tnemate: farsi i 


lagtonivg mre SITY ae aatoniogg ibis . a 


ath honored? 2aher9090. tem, eal soi fone 2 9 hs ) 
‘bag: dtbiw bred-tTst eosi3 16D | 


\ 


“sashats sig eames) ges 1s 


bsitans ont zotdngech Heo oth np ine” 


resis Som vt xis niside seonte ns 
aap: oT | 


me) eine Fas Bi a 


-Lnebyo ae zsniduordse svoito elieo - 


$ Ee 128 
ne ee Me sn ay 


86. 


BAND1 - Solves for and writes displacements using a Gaussian 
elimination technique for symmetric banded matrices. 

STRESS - Computes nodal and elemental stresses, element 
principal stresses and directions and characteristic 
nodal stresses. 

3.2 Definition of Variables 


The major arrays and variables used in the code are defined 


below. 
AP - applied load vector 
DJ - determinate of the Jacobian 
EBM - matrix [B] 
ECM - constitutive matrix [C] 
ER - Young's modulus in the radial direction 
ESM - element stress matrix [C][B] 


ESMRC - matrix [C][B]2mr 

ESTIF - element stiffness matrix [K]© 

ET - Young's modulus in the tangential direction 
EZ - Young's modulus in the axial direction 


ISODIM - 1 for isotropic material, 2 for anisotropic material 


KODE - index of displacement or load conditions 
M - element number 
MAT - material number 


MBAND=MM = - half-band width 
NEQ = NN- number of equations 
NCW - number of corners of an element where weight 


compensation is required. 
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NODE 


RO 
SIGRO 
SIGRZO 
SIGTO 
SIGZO 
SRCHAR 
STCHAR 
ST LF 
SZCHAR 
TUCHAR 
U 


V 
WP 
WT 


shearing stress t 


nodal point number 
element node number 
number of materials 
number of elements 
number of nodal points 
Poissons ratio wv 


re 


Poissons ratio v 
YZ 


Poissons ratio Voy 
radial coordinate of node 

radial coordinate in an element 
characteristic radius 

unit material weight 

radial stress 

rz 

tangential stress 

axial stress 

characteristic radial stress 
characteristic tangential stress 
stiffness matrix [K] 
characteristic axial stress 
characteristic shearing stress 
radial displacement or force 
axial displacement or force 
nodal weight compensation term 


one quarter of the element weight 


Sy: 
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Z - axial coordinate of node 


ZCHAR - characteristic radius 


3.3 Guide for Data Input 

Gide Input 

a) Heading Card (1 card) (FORMAT 18A4) 
b) Problem Card (1 card)(FORMAT 316) 


NUMNP NUMEL NUMAT 


c) Characteristics card (1 card)(FORMAT 6F12.4) 


RCHAR ZCHAR SZCHAR SRCHAR STCHAR TUCHAR 


d) Material cards (1 for each material) (FORMAT 3F12.0, 


- 4F6.0, 16) 
ER ET) EZ PRRT PRRZ PRTZ RO ISODIM 
e) Nodal point cards (FORMAT 216, 4F12.0) 
NODE KODE R v4 U V 


KODE 0 0 = R FORCE AND Z FORCE 
1 1 = R DISPLACEMENT AND Z DISPLACEMENT 
R DISPLACEMENT AND Z FORCE 


10 


O 1 = R FORCE AND Z DISPLACEMENT 
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RAHOUT = AAMT: 
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Nodal point cards must be in order. 


If nodes are omitted, the program generates co-ordinates for 


intermediate nodal points and sets 


KODE = 00, U = 0, V = O for these points. 


f) Element Cards (FORMAT 616) 


NP(1,M) 


NP(2,M) NP(3,M) NP(4,M) MAT(M) 


Element cards must be in order. If elements are omitted, the 


program will generate intermediate elements given a starting 


and ending element if node numbers along corresponding sides 


increase steadily by one. 


(ii) 


Progam Limitations 


Number of nodal points < 400 

Number of elements < 400 

Half-band width < 80 

Number of materials < 8 

Elements may not be used where their centroids lie on 


or outside their boundary. 
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s.4 Proqram Flow Chart 
MAIN 
\/ 
READIN 
\] 
— Data errors 
yes a 


| ELSTIe =a RETA 
y 


ASTIF no \/ 
— egative ele- 
MODIFY gent volume? 
yes\ 


> 
STRESS ee BETA 


V 


Go to 
next 
problem 


va 
WAY 


Is there a zero or negative entry on the main 
diagonal of the triangularized matrix? 


Figure 37. Program Flow Chart. 
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_ 3.5 PROGRAM 


COMMON ER(8),ET(8) ,EZ(8) ,PRRT(8) ,PRRZ(8) ,PRTZ(8), 
1 2(400),U (400) ,V (400), ISODIM (8) ,STIF(800,80) 
COMMON ECM (4,4),EBM(4, 8), ESM (4,8) ,WT,NUMNP,NUMEL, 
1 NP(4, 400) ,MAT (400) ,NEQ, NBAND, M, LM (8) 

COMMON RCHAR,ZCHAR,SZCHAR,SRCHAR,STCHAR,TUCHAR, 

1 RR(400) , ZR (400) ,SIGZO (400) ,SIGRO(400) ,SIGTO (400) 
COMMON RO(8),R(400),AP (800), ESTIF (8,8), 

INUMAT, KODE (400) , SIGRZO (400) 


Les 
1 CALL READIN 
Cc 

CALL ASTIF 
Cc 

CALL BAND1 
Cc 

CALL STRESS 
c 


GO TO 1 
END 
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SUBROUTINE READIN 


THIS SUBROUTINE READS AND PRINTS MATERIAL DATA, 
NODAL DATA, ELEEMENT DATA. IT GENERATES 
COORDINATES OF INTERMEDIATE NODAL POINTS AND 
CALCULATES THE BAND WIDTH AND NUMBER OF 
EQUATIONS 


COMMON ER(8),ET(8) ,£Z(8) ,PRRT(8) ,PRRZ(8) ,PRTZ(8), 
1 Z (400) ,U(400) ,¥ (400) , LSODIM(8) ,STIF(800,80) 
COMMON ECM(4,4),EBM(4,8),ESM(4,8),WT,NUMNP, NUMEL, 
1 NP (4,400), MAT (400) ,NEQ, MBAND,M, LM (8) 

COMMON RCHAR,ZCHAR,SZCHAR, SRCHAR,STCHAR, TUCHAR, 

1 RR(400) ,Z2R (400) ,SIGZO (400) , SIGRO(400) ,SIGTO (400) 
COMMON RO(8),R (400) ,AP (800) , ESTIF(8,8) , 

INUMAT, KODE (400) , SIGRZO (400) 


DIMENSION HED(18) 
READ PRELIMINARY INFORMATION 

READ(5,1000) HED,NUMNP,NUMEL,NUMAT 
WRITE(6,2000) HED, NUMNP,NUMEL, NUMAT 
READ (5,4000) RCHAR, ZCHAR,SZCHAR,S RCHAR, STCHAR, 
1TUCHAR 
WRITE(6,4010) RCHAR,ZCHAR,SZCHAR,SRCHAR,STCHAR, 
ITUCHAR 

4000 FORMAT (6F12.4) 

4010 FORMAT (//,*R CHARACTERISTIC =",F12.6/1H ,"Z% CHAR 
JACTERISTIC =",F12.6/1H ,*SIGMA Z CHARACTERISTIC 
2 =", F12.6/1H ,*SIGMA R CHARACTERISTIC =',F12.6/ 
31H ,'SIGMA T CHARACTERISTIC =',F12.6/1H , 
4 ‘SIGMA RZ CHARACTERISTIC =",F12.6// ) 


READ AND WRITE MATERIAL PROPERTIES 


WRITE (6, 2005) 
DO 10 M=1,NUMAT 
READ (5,1010) ER(M) ,ET(M),EZ(M),PRRT(M) ,PRRZ(M), 
1PRTZ(M) ,RO(M) , ISODIM(M) 

10 WRITE(6,2010) M,ER(M),ET(M),EZ(M),PRRT(M),PRRZ(M), 
IPRTZ (M) ,RO(M) , LSODIM {M) 


READ AND WRITE NODAL DATA AND GENERATE INTERMEDIATE 
NODAL DATA 


WRITE (6, 2014) 
WRITE (6, 2015) 
L=1 
READ (5,1020) N,KODE(N) ,R(N) ,Z(N) ,U(N),V(N) 
GO TO 40 

20 READ(5,1020) N,KODE(N) ,R(N) ,Z(N) ,U(N) /V(N) 
DN = N-L 
DR=(R(N)-R(L)) /DN 
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DZ= (Z(N)-Z(L)) /DN 
L=L+1 


IF(N-L) 50,40, 30 
30 R{(L)=R{L-1)+DR 

Z(L) =Z (1-1) +DZ 

KODE(L)= 0 

DEL) f= =O 

V(L)= 0 

GO TO 25 


40 WRITE(6,2020) N,KODE(N),R(N) ,Z(N),U(N) ,V(N) 
IP (NUMNE-N) 50,60, 20 


50 WRITE (6,2025) N 
CALL EXIT 


60 WRITE(6, 2016) 
WRITE (6, 2015) 


WRITE (6,2020) (N,KODE(N) ,R(N) ,Z(N) ,U(N) ,V(N), 


1N=1, NUMNE) 


READ AND WRITE ELEMENT DATA 


AAA 


WRITE (6, 2031) 
WRITE(6, 2030) 
NL=0 
51 IF (ML.GE.NUMEL) GO TO 70 


READ (5,1035) M,NP(1,M) ,NP(2,M) ,NP(3,M) ,NP(4,M), 


IMAT (M) 


WRITE (6,2035) M,NP(1,M) ,NP(2,M) ,NP(3,M) ,NE(4,M) , 


1MAT (M) 
MM=ML+1 
IF (MM.EQ.M) GO TO 65 


MM1=M-1 
DO 63 I=MM,MM1 
DO 62 J=1,4 
62 NP {J,I)=NP(J,I-1) +1 
MAT (I) =MAT {I-1) 
63 CONTINUE 


65 ML=M 
GO TO 51 
70 CONTINUE 
WRITE (6, 2032) 
WRITE (6, 2030) 
WRITE(6,2035) (M,(NP(J,M) p-J=1,4) ,MAT(M), 
1IM=1, NUMEL) 


GY@i-G? 


L=0 
DO 80 S8=1, NUMEL 


DETERMINE BAND WIDTH AND NUMBER OF EQUATIONS 


93% 
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DO 80 I=1,3 
II=1I+1 
DO 80 J=I1,4 
K= LABS (NP (I,M)-NP (J,M)) 
IF (K.GT.L) L=K 
80 CONTINUE 


MBAND = 2*(L¢1) 
NEQ= 2*NUMNP 


Cc 
WRITE (6, 2040) MBAND,NEQ 
IF (MBAND.LE. 80.AND.NEQ.LE.800) GO TO 90 
WRITE(6, 2050) 
CALL EXIT 

c 

90 WRITE (6, 3000) 

3000 FORMAT( * READIN COMPLETED * ///) 
RETURN 

¢ 


C FORMAT STATEMENTS 
1000 FORMAT (18A4/ 316) 
2000 FORMAT (1H1,10X,18A4,//// 


TAH aG 26H NUMBER OF NODAL POINTS = ,16/ 
2 RHE 26H NUMBER OF ELEMENTS = ,16/ 
SAAHES 26H NUMBER OF MATERIALS = ,I6) 


2005 FORMAT (///,1H ,10X, 21H MATERIAL PROPERTIES // 
11X,8SHMAT.NO. ,4X,13HMODULOUS E(R),6X,4HE(T) ,6X,4 
2HE (Z) ,4X,21HPCISSONS RATIO PR(RT),5X,6HPR(ZR), 
35X,6HPR(ZT) ,4X,11HUNIT WEIGHT,4X,6HISODIM ,//) 

1010 FORMAT (3F12.0,4F6.C,I6) 

2010 FORMAT (1H ,15,F20.0,2F10.0,F24.3,2F11. 3,F15.2,4X 
1,16) 

2014 FORMAT( '1*, 5X, "OUTPUT OF INPUT NODAL DATA *) 

2015 FORMAT (///,10X,19H NODAL POINT OUTPUT, /// 
11H ,* NODE KODE R CORD Z COORD R F 
2ORCE Z FORCE'//) 

2016 FORMAT (°1*,5X,*OUTPUT CF COMPLETE NODAL DATA") 

1020 FORMAT (216,4F12.0) 

2020 FORMAT (14,16, F13.6,3F12. 6) 

2025 FORMAT (1HO,28H ERROR IN NODAL DATA,NODE = ,I4) 

2030 FORMAT (///,10X, 13H ELEMENT DATA ///, 
1° ELEM I 3 K L MAT. NO.*//) 

2031 FORMAT (*1*,5X, "OUTPUT OF INPUT ELEMENT DATA* ) 

2032 FORMAT (*1",5X, "OUTEUT CF COMPLETE ELEMENT DATA‘) 
2035 FORMAT (I4, 516) 
1035 FORMAT (616) 
2040 FORMAT (///10X, 22H BAND WIDTH = 5 16/ 
1 10X, 22H NUMBER OF EQUATIONS = ,1I6) 
2050 FORMAT (///10X,*PROBLEM EXCEEDS SPECIFIED LIMITS’) 
G 
END 
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SUBROUTINE ASTIF 


THIS SUBROUTINE TAKES EACH ELEMENT IN TURN AND 
FORMS THE ELEMENT STIFFNESS MATRIX (BY CALLING 
ELSTIF). IT ASSEMBLES THE ELEMENT STIFFNESSES INTO 
STIF, ASSEMBLES THE APPLIED LOAD VECTOR (AP), AND 
MODIFIES THE ASSEMBLAGES FOR DISPLACEMENT 

BOUNDARY CONDITIONS (BY CALLING MODIFY) 


AAAAANANAANA 


COMMON ER(8),ET(8) ,EZ(8) ,PRRT(8) ,PRRZ(8) ,PRTZ(8) , 
1 2(400) ,U (400) ,¥ (400) , ISODIM (8) , STIF(800,80) 
COMMON ECM (4,4),EBM(4, 8), ESM (4,8) ,WT,NUMNP,NUMEL, 
1 NP (4,400) ,MAT (400) ,NEQ,MBAND,M, LM (8) 

COMMON RCHAR, ZCHAR,SZCHAR, SRCHAR,STCHAR,TUCHAR, 

1 BR(400) ,ZR (400) ,SIGZO (400) , SIGRO(400) , SIGTO (400) 
COMMON RO(8) ,R(400),AP (800), ESTIF(8,8) , 

INUMAT, KODE (400) , SIGRZO (400) 
COMMON/SLOP/WP (800) , NCW 


INITIALIZE APPLIED LOAD VECTOR AND MASTER 
STIFFNESS MATRIX AND ECM 


ORONO 


DO 10 I=1,NEQ 
WP (I) =0.0 
AP (I)=0.0 
DO 10 J=1,MBAND 
10 STIF (1,J)=0.0 
DOeZ207 I= 1, 4 
DO 21 J=1,8 
21 EBM (I,J) =0.0 
HO. 20 d= 174 
20 ECM (I, J) =0.0 
C FORM ELEMENT CONSTITUTIVE MATRIX (ECM) IF NUMAT=1) 


IF(NUMAT.NE.1) GO TO 30 
IF (ISODIM(1).NE.1) GO TO 22 
PR=PRRT (1) 
COM=ER (1) *(1.-PR) /((1-+PR) *(1.-2.*PR) ) 
COM1=PR/(1.-PR) 
ECM (1,1) =COM 
ECM (1, 2) =COM1*COM 
ECM (1, 3) =COM1*COM 
ECM (2, 1) =COM1*COM 
ECM (2,2) =COM 
ECM (2, 3) =COM1*COM 
ECM (3, 1) =COM1*COM 
ECM (3, 2) =COM1*COM 
ECM (3, 3) =COM 
ECM (4,4) =COM* (1.-2.*PR) /(2.*(1.-PR8) ) 
GO TO 30 
Cc ANISOTROPIC CASE 

22 EPRRT=PRRT (1) 
EPRZT= PRTZ (1) 
EPRZR=PRRZ (1) 
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VET=ET (1) 

VER=ER (1) 

VEZ=EZ (1) 

ERZ=VER/VEZ 

ETR=VET/VER 

ETZ=VET/VEZ 

COM=1./ (1. ~EPRRT ** 2* FT R- EPRZR*¥* 2* ER Z-EPRZT 2 
1* ETZ-2.* EPRZR*EPRZT*EPRRT*ETZ) 
G=VEZ/2.¥*(1. +EPRZR) 

ECM (1, 1) =VEZ* (1. -EPRRT**2*ETR) *COM 

ECM (1, 2) = (EPRZR*VER+EPRRT*EPRZT* VET) *COM 
ECM (1, 3) =VET* (EPRZR*EPRRT+EPRZT) *COM 

ECM (2, 2) =VER* (1.-EPRZT**2*ETZ) *COM 

ECM (2, 3) =VET* (EPRRT+ EP RZR*EPRZT*ERZ) *COM 
ECM (3, 3)=VET* (1. -EPRZR**2*ERZ) *¥COM 
ECM (4, 4) =G 

ECM (2, 1) =ECM (1,2) 

ECM (3, 1) =ECM (1, 3) 

ECM (3, 2) =ECM (2, 3) 


30 DO 45 s=1,NUMEL 
CALL ELSTIF 
c 
Cc ASSEMBLE ELSTIF INTO MASTER STIFFNESS MATRIX 
Cc 


DO 35 I=1,4 

12=2*1 

LM (12) =2*NP (I, M) 
35 LM (12-1) =LM (12) -1 


DO 40 I=1,8 

II=LM (1) 

I1=14+3¥* (1/2) ¢4* (-1/7+1/6-1/5-I/3) 

DO 40 J=1,8 

JJ=LM (J) -1I+1 

J1=39 +3 * (J/2) +4*(-3/74+I/6-J/5-J/3) 

IF (JJ.LE.0) GO TO 40 

STIF (11,JJ) =STIF (11,JJ) +ESTIF (11,J1) 
0 CONTINUE 


4 
Cc 
G ADD GRAVITY LOADS INTO WP VECTOR 
C 


DO 45 I=2,8,2 
II=LM (I) 
WP (II) =WP (IL) -wT 
45 CONTINUE 
fe 
C MULTIPLY NODAL LOADS BY 2*PI*R AND 
C ADD NODAL LOADS & GRAVITY LOADS INTO AP VECTOR 
Cc 


DO 50 N=1,NUMNP 
IF (KODE(N) -1) 46,47,48 
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V (N)=V (N) ®2.*3.1416*ABS (R(N) ) 

U (N) =U (N) #2. *3.1416*ABS (R(N) ) 
48 CONTINUE 

IF (KODE(N)-10) 51,49,51 

V (N)=V (N) #2. *3.1416¥*ABS (2 (N) ) 
51 N2=2*N 

AP (N2) =WP(N2) +V¥(N) 

AP (N2-1) =AP (N2-1) +U (N) 


MODIFY STIFFNESS AND LOAD VECTOR FOR DISPLACEMENT 
BCUNDARY CONDITIONS 


NCW=0 

DO 102 N=1,NUMNP 

N2=2*N 

II=N2-1 

IF (KODE (N)-1) 100, 80,60 
60 CALL MODIFY (II,N) 

IF (KODE(N).NE.11) GO TO 100 
80 CALL MODIFY (N2,N) 

NC W= NC W411 

GO TO 102 
100 WP (N2) =0.0 
102 CONTINUE 


RETURN 
END 
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SUBROUTINE ELSTIF 


C THIS SUBROUTINE FORMS THE ELEMENT STIFFNESS MATRIX 
(ESTIF) FCR LINEAR STRAIN TRAPEZOID 


C 


c 


COMMON ER(8),ET(8) ,EZ(8),PRRT(8) ,PRRZ(8) ,PRTZ(8), 
1 Z(400),U (400) ,V¥ (400) , ISODIM(8) , STIF(800,80) 
COMMON ECM(4,4),EBM(4,8) ,ESM (4,8) ,WT,NUMNP,NUMEL, 
1 NP (4,400) , MAT (400) ,NEQ, MBAND, M, LM (8) 

COMMON RCHAR,ZCHAR,SZCHAR,SRCHAR,STCHAR,TUCHAR, 

1 RR(400) ,2R(400) ,SIGZO (400) , SIGRO(400) , SIGTO (400) 
COMMON RO(8),R (400) ,AP (800) , ESTIF(8,8) , 

INUMAT, KODE (400) ,SIGRZO (400) 

COMMON/BETTER/ R&T (4) ,Z2T(4),FIJ (4,4) ,DJ,RAD, 
1ESMRC (4, 8) 


GP = 1.0/SQRT (3.0) 


ACCUMULATES ELEMENT COORDINATE LOCATIONS 


20 


DO 20 °I=1,4 
J=NP (I,M) 
RT (I)=R(J) 
ZT (I)=Z(d) 


FORM CONSTITUTIVE MATRIX 


Ze 


NM=1 
IF (NUMAT.EQ.1) GO. TO 30 
NM=MAT (M) 
IF {ISODIM(NM).NE.1) GO TO 22 
PR=PRRT (NM) ) 
COM=ER (NM) *(1.-PR) /{(1.+PR) * (1.-2. *PR) ) 
COM1=PR/ (1.-PR) 
ECM (1,1) =COM 
ECM (1, 2) =COM1*COM 
ECM (1,3) =COM1*COM 
ECM (2, 1) =COM1*COM 
ECM (2,2) =COM 
ECM (2, 3) =COM1*COM 
ECM (3, 1)=COM1*COM 
ECM (3, 2)=COM1*COM 
ECM (3, 3) =COM 
ECM (4, 4) =COM* (1.-2.*PR) / (2. *(1.-PR) ) 
GO TO 30 
ANISOTROPIC CASE 
EPRRT=PRRT (NM) 
EPRZT= PRTZ (NM) 
EPRZR=PRRZ (Ni) 
VET=ET (NM) 
VER=ER (NM) 
VEZ=EZ (NM) 
ERZ=VER/VEZ 
ETR=VET/VER 
ETZ=VET/VEZ 
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COM=1. / (1. -EPRRT** 2* FT R- EPRZR**2*ERZ- EPRZT **2 
1*ETZ—2.* EPRZR¥EPRZT*EPRRT*ETZ) 
G=VEZ/2.*(1.+EPRZR) 
ECM (1, 1) =VEZ* (1. -EBRRT**2*ETR) *COM 
ECM (1, 2) = (EPRZR*VER*+EPRRT*EPRZT* VET) *COM 
ECM (1, 3) =VET* (EPRZR*¥EPRRT+EPRZT) *COM 
ECM (2, 2) =VER* (1.-EPRZT**2*ETZ) *COM 
ECM (2, 3) =VET* (EPRRT#+EPRZR*EPRZT*ERZ) *COM 
ECM (3, 3) =VET* (1.-EPRZR**2*ERZ) *COM 
ECM (4,4) =G 
ECM (2, 1)=ECM (1,2) 
ECM (3, 1) =ECM (1, 3) 
ECM (3, 2) =ECM (2,3) 
ZERO'S STIFFNESS MATRIX 
DO 19 I=1,8 
pO 19 J=1,8 
ESTIF (1,3)=0.0 
VOL=0. 0 


FORM ELEMENT STIFFNESS MATRIX ({ESTIF) 


DOES INTEGRATION QUADRANT BY QUADRANT 
USING GAUSSE LEGENDRE QUADRATURE 


28 


DO 28 K=1,4 
K1=- 14 2* (K/2-2* (K/4) ) 
K2=- 1+ 2* (K/3) 

EPS= FLOAT (K1) *GP 
PNU= FLOAT (K2)¥*GP 


CALL BETA(EPS, PNU) 
VOL=VOL+DJ 

DO 28 I=1,8 

DO 28 J=1,8 

DO 28 L=1,4 

ESTIF (I,J) =ESTIF (I,J) #EBM (L,I) *ESMRC (L,J) 
WT=VOL*3.1416*RAD*RO (NM) /2.0 

RETURN 


END 
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SUBROUTINE BETA (EPS, PND) 

COMMON ER(8),ET(8) ,EZ(8),PRRT(8) ,PRRZ(8) ,PRTZ (8), 
1 Z(400),U (400) ,V (400) , ISODIM (8) , STIF(800,80) 
COMMON ECM (4,4) ,EBM(4, 8), ESM (4,8) ,WT,NUMNP, NUMEL, 
1 NP (4,400) , MAT (400) ,NEQ, MBAND, M, LM (8) 

COMMON RCHAR,ZCHAR,SZCHAR, SRCHAR,STCHAR,TUCHAR, 
1 RR(400) ,2ZR (400) ,SIGZO (400) ,SIGRO(400) , SIGTO (400) 
COMMON RO(8),R (400) ,AP (800) , ESTIF(8,8) , 

1NUMAT, KODE (400) , SIGRZO (400) 

COMMON/BETTER/ RT(4),2T(4),FIJS (4,4) ,DJ,RAD, 
1ESMRC (4, 8) 

R42=RT (4) -RT (2) 

Z242=Z2T (4) -ZT (2) 

R31=RT (3) -RT (1) 

Z31=2ZT (3) -ZT(1) 

RE34= (RT (3) -RT (4)) *EPS 

ZE34= (ZT (3) -ZT (4)) *EPS 

RE12= (RT (1) -RT(2)) *EPS 

ZE12=(ZT(1)-ZT (2)) *EPRS 

RN4&1= (RT (4)-RT(1)) *PNU 

ZN41= (ZT (4) -ZT (1) ) *PNU 

RN23= (RT (2) -RT(3)) *PNU 

ZN23= (ZT (2) -ZT (3) ) *PNU 

EBM (1, 5) =(R42+RE34+4RN23) /8.0 

EBM (1,6) = (RN4 1-R31-RE34) /8.0 

EBM (1, 7) =(R42+RE12+BN41) / (-8.0) 

EBM (1, 8) = (R31+RE12-RN23) /8.0 

EBM (2, 1) =- (Z42+ZE34+2ZN23) /8.0 

EBM (2, 2)=(Z31+ZE34-ZN41) /8.0 

EBM (2, 3) =(Z42+ZE12+ZN41) /8.0 

EBM (2, 4) = (ZN23-231-ZE12) /8.0 

Dd=6. 


FORM ELEMENT B MATRIX (EBM) 


POl 2Se1=4,.48 
DJ=DJ+ EBM (1,1+#4) *ZT (I) 

EBM (4, I) =EBM(1,1+4) 

25 EBM(4,1I+4) =EBM (2,1) 

if 20d 20E2 00) GO-20 75 
RAD=.25*((1.-EPS) * (1.-PNU) #RT (1) +(1.+EBPS) * 
1(1.-PNU) ¥RT (2) #(1. +EPS) *(1.+PNU) *RT (3) + 
2(1.-EPS)*(1.#PNU) *RT (4) ) 

RAD=ABS (RAD) 

EBM (3, 1) =Dd* (2 25*(1.-EPS) *(12-PNU)) /RAD 
EBM (3,2) =DJ* (.25*(1.+EPS) *(1.-PNU)) /RAD 
EBN (3,3) =DJ¥*(.25*(1.+EPS) *(1.#PNU) ) /RAD 
EBM (3, 4) =Dd* (. 25*(1.-EPS) *(1.#PNU) ) /RAD 


FORM ELEMENT STRESS MATRIX (ESM) 
DO 27 I=1,4 


DO QV7id=t)3 
ESM(I,J)=0.0 
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DO 26 L=1,4 
26 ESM(I,J)=ESM(I,J)+ECM(I,L) *EBM(L,J) 
ESM (I, J) =ESM (I,J) /DJ 
27 ESMRC (I,J) =ESM (I,J) *2. *3. 1416*RAD 
RETURN 
75 WRITE(6,1000) M 
1000 FORMAT (1H1,*VOLUME OF ELEMENT*,14,°IS LE 
1SS THAN ZERO‘) 
CALL EXIT 
80 RETURN 
END 
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SUBROUTINE BAND1 


THIS SUBROUTINE SOLVES FOR DISPLACEMENTS USING A 
GAUSSIAN ELIMINATION TECHNIQUE FOR SYMMETRIC 
BANDED MATRICES STOREC IN CORE 


COMMON ER(8),ET(8) ,EZ(8),PRRT(8) ,PRRZ(8) ,PRTZ(8), 
1 Z(400),U (400) ,¥(400),ISODIM(8), A(800,80) 
COMMON ECM(4,4),EBM(4,8),ESM(4,8),WT,NUMNP,NUMEL, 
1 NP (4,400) ,MAT (400), NN,MM,M,LM(8) 

COMMON RCHAR,ZCHAR,SZCHAR,SRCHAR,STCHAR,TUCHAR, 

1 RR(400) ,ZR (400) ,SIGZO (400) , SIGRO(400) , SIGTO (400) 
COMMON RO(8),R (400), B(800),ESTIF(8,8) , 

INUMAT, KODE (400) , SIGRZO (400) 

DIMENSION NODNO (800) 


TRIANGULARIZE AND REDUCE RIGHT HAND SIDE 
NL=NN-MM+1 
NM=NN-1 
MR=MM 


DO 100 N=1,NM 
IF (A(N,1).LE.0.) GO TO 700 
BN=B (N) 

B(N)=BN/A (N, 1) 

IF (NeGT.NL) MR=NN-N#1 


DO 100 L=2,MR 
IF (A(N,L).~EQ.0.)GO TO 100 
= A(N,L)/A(N,1) 


= N¢L-1 

= 0 

DO 50 K=L,¥M@R 
J=J+1 


A{I,J)= A(1I,J) -C*¥A(N,K) 
B (1) =B {I) -C*BN 
A(N,L) =C 

0 CONTINUE 


BACK SUBSTITUTE 
I=NN 


B(NN)=B(NN) /A(NN,1) 

DO 600 N=1,NM 

I=I-1 

IF (N.LT.MM) MR= N¢+1 
DO 600 J=2,MR 

K=I+J-1 


600  B(I) =B{I)- A(I,d) *B(K) 


DO 650 N=1,NN 
650 NODNO(N)=1+4 (N-1) /2 
WRITE(6,2001) (NODNO(I),% (I) -I=1,NN) 


2004 ‘FORMAT (21%, , * NODAL DISPLACEMENTS * 47/7 
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18 NODE*,11X,°U*,6X,"NODE",11X,"VI// 
22 (16.,815~.5)} 
RETURN 

700 WRITE (6,2000) N 
GRELYEXIT 

2000 FORMAT (1HO,*ZERO OR NEGATIVE ELEMENT ON MAIN DIA 

1GONAL OF TRIANGULARIZED MATRIX FOR EQUA 
2TION?, 15) 


END 
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SUBROUTINE MODIFY ({1I,N) 


THIS SUBROUTINE MODIFIES KSTIF AND AP IF A 
DISPLACEMENT BOUNDARY CONDITION IS IMPOSED IN 
EQUATION I ASSOCIATED WITH NODAL POINT N 


COMMON ER(8),ET(8) ,EZ(8),PRRT(8) ,PRRZ(8) ,PRTZ(8), 
1 Z (400) ,U (400) ,V (400) , LSODIM (8) , STIF (800,80) 
COMMON ECM (4,4),EBM(4, 8), ESM (4,8) ,WT,NUMNP,NUMEL, 
1 NP(4, 400) , MAT (400) ,NEQ,MBAND,M, LM (8) 

COMMON ECHAR,ZCHAR,SZCHAR,SRCHAR,STCHAR,TUCHAR, 

1 RR(400) , ZR (400) ,SIGZO (400) ,SIGRO(400) ,SIGTO (400) 
COMMON RO(8),R(400),AP (800), ESTIF(8,8) , 

INUMAT, KODE (400) , SIGRZO (400) 


DISP=U (N) 
IF ((I-2*N).EQ.0) DISP=V(N) 


DO 50 J=2, MBAND 

IL=I+J-1 

IU=I-J+1 

IF{IU.LE.0)GO TO 10 

AP (10) =AP(IU) - STIF(IU,J) *DISP 
STIF(IU ,J)=0.0 

IF(IL.GT.NEC) GO TO 50 

AP (IL)= AP({IL)- STIF (I,J) *DISP 
STIF (I,J) =0.0 

CONTINUE 

AP(I)=DISP 

STIF (I, 1)=1.0 

RETURN 

END 


108. 


SUBROUTINE STRESS 


THIS SUBROUTINE FORMS THE LEMENT STRESS MATRIX 
(ESM), MULTIPLIES BY THE FLEMENT DISPLACEMENT 
VECTOR (ELDISP) AND RECORDS THE STRESSES IN 
SIGEL. IT THEN COMPUTES THE PRINCIPAL STRESSES 
AND DIRECTIONS (SIGP) 


ARQANANN 


COMMON ER(8),ET(8) ,EZ(8) ,PRRT(8) ,PRRZ(8) ,PRTZ(8), 
1 Z(400) ,U(400) ,V (400), ISODIM (8) ,STIF (800,80) 
COMMON FCM (4,4),EBM(4,8) ,ESM (4,8) ,WT,NUMNP,NUMEL, 
1 NP(4, 400) , MAT (400) ,NEQ, MBAND,M, LM (8) 

COMMON KCHAR,ZCHAR,SZCHAR, SRCHAR, STCHAR, TUCHAR, 

1 RR(400) ,ZR (400) , SIGZO (400) , SIGRO(400) ,SIGTO (400) 
COMMON RO(8),R(400),AP (800), ESTIF(8,8), 

INUMAT, KODE (400) , SIGRZO (400) 


DIMENSION SIGEL(400,4) ,SIGP(400,7) , ELDISP(8), 
1SIGA(400,4) ,KOUNT (400) 

COMMON/BETTER/ RT(4),ZT(4),FId (4,4) ,DJ,RAD, 
1ESMRC (4, 8) 

COMMON/SLCE/WP (800) , NCW 


ZERO'S SIGEL,KOUNT & SIGA 


Cy ALea 


DO 1 I=1,NUMEL 
DO 1 J=1,4 
1 SIGEL{I,J)=0. 
DO 5 N=1,NUMNP 
KOUNT(N) =0 
DO 5 J=1,4 
5 SIGA (N,J)=0.0 


DO 15 M=1,NUMEL 
COMPUTE ELEMENT DISPLACEMENTS 


aA 


poveosieyys4 
J=NP (I ,M) 
ELDISP {1) =AP (2*J-1) 
ELDISP (1+4) =AP (2*J) 
KOUNT (J) =KOUNT (J) +1 
RT (I)=R (J) 
10 2948) =2 (3) 
DO 15 K=1,4 
L=NP (K , M4) 
EPS=FLCAT (-14+2* (K/2-2%* (K/4)) ) 
PNU=FLOAT (- 1+2* (K/3) ) 
CALL BFTA(EEPS, PNU) 
DOPAS 1=1,4 
Cc 
C ACCUMULATE FOR NODAL STRESSES 
c 
DO 15 J=1,8 
15 SIGA (L,I) =SIGA (L,I) tESM(1,J) *ELDISP (J) 
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Cc 
C FIND AVERAGE NODAL STRESSES 
Cc 

DO 110 K=1,NUMNP 

RK=KOUNT (5) 

DO 110 I=1,4 
110 SIGA(N,1I)= SIGA(N,I) /RK 
‘o 
C COMPUTE ELEMENT STRESSES 

DO 16 #=1,NUMEL 

DO 16 K=1,4 

L=NP {K,) 

DO 16 I=1,4 

16 SIGEL(#,1)=SIGEL (4,1) #SIGA({L,I) /4. 

Cc 
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C COMPUTE ELEMENT PRINCIPAL STRESSES AND DIRECTIONS 


Cc 
DO 100 s=1,NUMEL 
SIGM=(SIGEL(M,1) +SIGEL (M,2)) /2. 
SIGD2= (SIGEL (M, 1) -SIGEL(M, 2) ) /2. 
RAD = SQRT(SIGD2**2 #SIGEL(,4%) **2) 
SIGP (N,1)=SIGMN#RAD 
SIGP (§,2)= SIGH-RAD 
SIGP (4,3) =SIGEL {*, 3) 

SIGP (M,4)= (SIGP (EH, 2) -SIGP (M, 3)) /2. 
SIGP (N,5)=(SIGP (#8, 1) -SIGP? (M,3))/2. 
SIGP (8 ,6)=(SIGP(M, 1) -SIGP(M,2))/2. 


SIGP (8,7) =0.5*57.29578*ATAN2 (SIGEL(N,4%) ,SIGD2) 


100 CONTINUE 


FIND MAXINUN ELEMENT STRESSES 


AAA 


S1IG1=SIGP (1,1) 
SIG2=SIGP (1,2) 
SIG3=SIGP (1,3) 
SIG&=SIGP (1,4) 
SIG5=SIGP {1,5) 
SIG6=SIGP (1,6) 
™1=0 
M2=0 
“3=0 
Mu=0 
M5=0 
M6=0 
DO 145 MN=1,NUMEL 
IF (SIGP(N,1).LT.S1G1) GO TO 115 
SIG1=SIGP (M, 1) 
M1=48 
115 IF (SIGP(M,2).GT.SIG2) GO TO 120 
SIG2=SIGP (M, 2) 
42= 
120 IF {(SIGP(M,3).LT.SIG3) GO TO 130 
SIG3=SIGP (é, 3) 
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M3=M 

130 IF (SIGP(M,4).LT.~SIG4) GO TO 135 
SIG4=SIGE (M,4) 
M4=M 

135 IF (SIGP(M,5).LT.SIG5) GO TO 141 
SIG5=SIGP (M,5) 
MS=M 

141 IF (SIGP(M,6).LT.SIG6) GO TO 145 
SIG6=SIGP (M,6) 
M6O=M 

145 CONTINUE 
C FIND CHARACTERISTIC NCDAL VALUES 

DO 300 N=1,NUMNP 
RR (N)=R(N) /RCHAR 
ZR (N)=Z(N) /ZCHAR 
SIGZO(N) =SIGA(N, 1) /SZCHAR 
SIGRO(N)=SIGA(N, 2) /SRCHAK 
SIGTO(N) =SIGA (N, 3) /STCHAR 
SIGRZO (N) =SIGA(N,4) /TUCHAR 

300 CONTINUE 
WRITE(6,2000) 
WRITE(6,2010) (M,(SIGEL(M,1) ,1=1,4),M=1,NUMEL) 
WRITE (6,2020) 
WRITE(6,2030) (M,(SIGP (M,I) ,I=1,7) »M=1,NUMEL) 
WRITE (6,2040)SIG1,M1,S1G2,M2,S1G3,M3,SIG4, M4, 
1SIG5,M5,SIG6,M6 
WRITE (6,2050) 
WRITE(6,2010) (N,(SIGA (N,1I) ,1=1,4) ,N=1,NUMNP) 
WRITE(6,3010) 
WRITE (6, 3020) (N,RR(N) ZH (N) »STGZO(N) ,SIGRO(N) , 
1SIGTO(N) ,SIGRZO(N) ,N=1,NUMNP) 

3010 FORMAT (1H1, "AVERAGE CHARACTERISTIC NODAL VALUES® 


V///1X, "NODE R CH. Z CH. SIGMA ZC 
2H. SIGMA R CH. SIGMA T CH. SIGMA R 
3Z CH.*//) | 


3020 FORMAT (I4,6F12.6) 

2000 FORMAT (1H1,10X, 21H R-Z ELEMENT STRESSES /// 
11X,°ELEM SIGMA Z SIGMA R SIGMA T 
2 SIGMA RZ'//) 

2010 FORMAT (I4,F11.6,3F 12.6) 

2020 FORMAT (1H1,10X,*ELEMENT PRINCIPAL STRESSES!/// 


11X,*° ELEM SIGMA I SIGMA Ti SIGMA If 
2 SIGMA IItl TAU RZ TAU RT TAU TZ 
3 SIGMA I DEG!//) 


2030 FORMAT (I4,F11.6,5F 12.6, F14. 3) 

2040 FORMAT (1H1,"MAXIMUM PRINCIPAL STRESSES IN RZ PL 
JANE= *,F11.6,' AND OCCURS IN ELEM*,I6//* MINIMU 
1M PRINCIPAL STRESS IN RZ PL 
2ANE= ',F11.6,* AND OCCURS IN ELEM*,I6//" MAXIMU 
3M HOOP STRESS =",F11.6,' AND OCCURS IN ELEM',I6/ 
G/* MAXIMUM TAU RZ =",F11.6,' AND OCCURS IN ELEM! 
5,16//" MAXIMUM TAU RT =*,F11.6," AND OCCURS IN E 
6LEM*,I6//" MAXIMUM TAU TZ =",F11.6," AND OCCUR 
7S IN ELEM',I6//) 
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108. 


2050 FORMAT (1H1, ‘AVERAGE NODAL STRESSES ee 
1* NODE SIGMA Z SIGMA R SIGMA a 


2GMA RZ'"//) 
WRITE (6,2060) 

2060 FORMAT (*0%,*REACTICN CORRECTION AT NODES DUE T 
10 WEIGHT‘) 
IF (NCW.EQ.0) GO TC 150 
WRITE (6,2070) 

2070 FORMAT (T10,*NODE* ,T20,*CORRECTION! ,//) 
T=0 
DO 140 N=1,NUMNP 
NN=2*N 
IF (WP (NN) ~EQ.0.0) GO TO 140 
WRITE (6,2080)N,WP (NN) 

2080 FORMAT (111,13,723,F10.4) 
I=1I+1 
IF (I1.GE.NCW) RETURN 

140 CONTINUE 
150 WRITE (6,2090) 

2090 FORMAT ('NONE') 
RETURN 
END 
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